/*----------------------------------------------------------------------*/
/* PROGRAM: treatment_effects.ado					*/
/*									*/
/* PURPOSE:								*/
/* [*] 	This code calculates the treated outcomes, untreated outcomes, 	*/
/*	and treatment effects using the baseline treated (BT) weights,  */
/*	the baseline untreated (BU) weights, the intervention treated 	*/
/*	(IT) weights, the intervention untreated (IU) weights, the 	*/
/*	randomized intervention sample treated (RIST) weights, the	*/
/*	randomized intervention sample untreated (RISU) weights, the 	*/
/*	local average (LA) weights, and the average (A) weights. The 	*/
/*	code below lists step-by-step how to derive these weights and 	*/
/*	the treatment effects using the weights. A good resource on 	*/
/*	what this code does is Table 1 of the paper.			*/
/*									*/
/* INPUTS:								*/
/* [*] 	arg_MTE_matrix:	The matrix containing the polynomial expression */
/*	for the MTE. This matrix should have one row and M+1 columns, 	*/
/*	where M is the order of the polynomial. The first term in the 	*/
/*	matrix is the intercept, the second the linear term, etc. The	*/
/*	code will work on any reasonable polynomial order M. This code 	*/
/*	can be applied to any MTE functional form that the user would 	*/
/*	like to use to derive the treatment effects (e.g., it can be	*/
/*	MTE(p), SMTE(p), or mte(p))					*/
/* [*] 	arg_MTO_matrix:	The matrix containing the polynomial expression */
/*	for the MTO. This matrix should have one row and M+1 columns, 	*/
/*	where M is the order of the polynomial. The first term in the 	*/
/*	matrix is the intercept, the second the linear term, etc. The	*/
/*	code will work on any reasonable polynomial order M. This code 	*/
/*	can be applied to any MTO functional form that the user would 	*/
/*	like to use to derive the treatment effects (e.g., it can be 	*/
/*	MTO(p), SMTO(p), or mto(p))					*/
/* [*] 	arg_MUO_matrix:	The matrix containing the polynomial expression */
/*	for the MUO. This matrix should have one row and M+1 columns, 	*/
/*	where M is the order of the polynomial. The first term in the 	*/
/*	matrix is the intercept, the second the linear term, etc. The	*/
/*	code will work on any reasonable polynomial order M. This code 	*/
/*	can be applied to any MTO functional form that the user would 	*/
/*	like to use to derive the treatment effects (e.g., it can be	*/
/*	MUO(p), SMUO(p), or muo(p))					*/
/* [*] 	arg_pB:	The macro variable containing the value of pB=P(D=1|Z=0)*/
/* [*] 	arg_pI:	The macro variable containing the value of pI=P(D=1|Z=1)*/
/* [*] 	arg_s_pI:The macro variable containing the value of s(pI)=P(Z=1)*/
/*									*/
/* OUTPUTS:								*/
/* [*] sc_BTTO: The treated outcome for baseline treated (scalar)	*/
/* [*] sc_BTUO:	The untreated outcome for baseline treated (scalar)	*/
/* [*] sc_BTTE:	The treatment effect for baseline treated (scalar)	*/
/* [*] sc_BUTO:	The treated outcome for baseline untreated (scalar)	*/
/* [*] sc_BUUO:	The untreated outcome for baseline untreated (scalar)	*/
/* [*] sc_BUTE:	The treatment effect for baseline untreated (scalar)	*/
/* [*] sc_ITTO:	The treated outcome for intervention treated (scalar)	*/
/* [*] sc_ITUO:	The untreated outcome for intervention treated (scalar)	*/
/* [*] sc_ITTE:	The treatment effect for intervention treated (scalar)	*/
/* [*] sc_IUTO:	The treated outcome for intervention untreated (scalar)	*/
/* [*] sc_IUUO:	The untreated outcome for intervention untreated(scalar)*/
/* [*] sc_IUTE:	The treatment effect for intervention untreated (scalar)*/
/* [*] sc_LATO:	The local average treated outcome (a scalar)		*/
/* [*] sc_LAUO: The local average untreated outcome (a scalar)		*/
/* [*] sc_LATE:	The local average treatment effect (a scalar)		*/
/* [*] sc_ATO:	The average treated outcome (a scalar)			*/
/* [*] sc_AUO:	The average untreated outcome (a scalar)		*/
/* [*] sc_ATE:	The average treatment effect (a scalar)			*/
/* [*] sc_RISTTO: The RIST treated outcome (a scalar)			*/
/* [*] sc_RISTUO: The RIST untreated outcome (a scalar)			*/
/* [*] sc_RISTTE: The RIST treatment effect (a scalar)			*/
/* [*] sc_RISUTO: The RISU treated outcome (a scalar)			*/
/* [*] sc_RISUUO: The RISU untreated outcome (a scalar)			*/
/* [*] sc_RISUTE: The RISU treatment effect (a scalar)			*/
/*									*/
/* EXAMPLE CALL TO PROCEDURE:						*/
/* matrix MTO = J(1,2,.)						*/
/* matrix MUO = J(1,2,.)						*/
/* matrix MTE = J(1,2,.)						*/
/*									*/
/* matrix MTO[1,1] = 10582						*/
/* matrix MTO[1,2] = -23601						*/
/* matrix MUO[1,1] = 3905						*/
/* matrix MUO[1,2] = -1383						*/
/*									*/
/* matrix MTE[1,1] = MTO[1,1] - MUO[1,1]				*/
/* matrix MTE[1,2] = MTO[1,2] - MUO[1,2]				*/
/*									*/
/* local pB = 0.152							*/
/* local pI = 0.411							*/
/* local s_pI = 0.344							*/
/*									*/
/* treatment_effects MTE MTO MUO  `pB' `pI' `s_pI'			*/
/*----------------------------------------------------------------------*/

capture program drop treatment_effects

program treatment_effects
	args arg_MTE_matrix arg_MTO_matrix arg_MUO_matrix arg_pB arg_pI arg_s_pI
	marksample touse
	
	* Print out diagnostics
	
	di in gr "Beginning calculation of treatment effects."
	
	if `arg_pB'==0 {
		di "WARNING: pB is equal to 0."
	}
	
	if `arg_pB'==1 {
		di "WARNING: pB is equal to 1."
	}
	
	if float(`arg_pB') == float(`arg_pI') {
		di "WARNING: pB is equal to pI."
	}
	
	if `arg_pI'==0 {
		di "WARNING: pI is equal to 0."
	}
	
	if `arg_pI'==1 {
		di "WARNING: pI is equal to 1."
	}
		
	* Calculate the uniform treatment effects in Mata
	
	local arg_s_pB = 1-`arg_s_pI'
	
	mata: calc_treat_eff("`arg_MTE_matrix'", "`arg_MTO_matrix'", "`arg_MUO_matrix'", `arg_pB', `arg_pI', `arg_s_pI', `arg_s_pB')		
end



/* Define the function for deriving the treatment effects in Mata*/
mata:
void calc_treat_eff(MTE_name, MTO_name, MUO_name, pB_name, pI_name, s_pI_name, s_pB_name)
{
	MTE_m = st_matrix(MTE_name)
	MTO_m = st_matrix(MTO_name)
	MUO_m = st_matrix(MUO_name)
		
	/* BT weights */

	/* Define the weights as 1 from 0 to pB */
	/* There is no reference to 0 and pB here, but later when we take the */
	/* sum of the weights we take it from 0 to pB */
	BT_w = 1
		
	/* Integrate the weights from 0 to pB */
	/* This is equivalent to calculating the sum of the BT weights from */
	/* 0 to 1 */
	BT_sum = polyeval(polyinteg(BT_w,1),pB_name) - polyeval(polyinteg(BT_w,1), 0)
		
	/* Normalize the weights to sum to 1 by dividing the weights by the */
	/* sum from 0 to 1 */
	polydiv(BT_w, BT_sum, BT_w_new, temp)
		
	/* Multiply MTO, MUO, and MTE by the weights calculated above*/
	w_BT_MTO = polymult(MTO_m, BT_w_new)
	w_BT_MUO = polymult(MUO_m, BT_w_new)
	w_BT_MTE = polymult(MTE_m, BT_w_new)
		
	/* Integrate MTO, MUO, and MTE from 0 to pB */
	/* If pB = 0, then evaluate the appropriate function at 0 and use that*/
	/* as the output */
	
	if (pB_name==0) {
		BTTO = polyeval(MTO_m,pB_name)
		BTUO = polyeval(MUO_m,pB_name)
		BTTE = polyeval(MTE_m,pB_name)	
	} 
	else {
		BTTO = polyeval(polyinteg(w_BT_MTO,1),pB_name) - polyeval(polyinteg(w_BT_MTO,1),0)
		BTUO = polyeval(polyinteg(w_BT_MUO,1),pB_name) - polyeval(polyinteg(w_BT_MUO,1),0)
		BTTE = polyeval(polyinteg(w_BT_MTE,1),pB_name) - polyeval(polyinteg(w_BT_MTE,1),0)
	}
			
	st_numscalar("sc_BTTO", BTTO)
	st_numscalar("sc_BTUO", BTUO)
	st_numscalar("sc_BTTE", BTTE)	
				
	/* BU weights */

	/* Define the weights as 1 from pB to 1*/
	/* There is no reference to 1 and pB here, but later when we take the */
	/* sum of the weights we take it from pB to 1 */
	BU_w = 1
		
	/* Integrate the weights from pB to 1 */
	/* This is equivalent to calculating the sum of the BU weights from */
	/* 0 to 1 */
	BU_sum = polyeval(polyinteg(BU_w,1),1) - polyeval(polyinteg(BU_w,1), pB_name)
		
	/* Normalize the weights to sum to 1 by dividing the weights by the */
	/* sum from 0 to 1 */
	polydiv(BU_w, BU_sum, BU_w_new, temp)
		
	/* Multiply MTO, MUO, and MTE by the weights calculated above*/
	w_BU_MTO = polymult(MTO_m, BU_w_new)
	w_BU_MUO = polymult(MUO_m, BU_w_new)
	w_BU_MTE = polymult(MTE_m, BU_w_new)
		
	/* Integrate MTO, MUO, and MTE from pB to 1*/
	/* If pB = 1, then evaluate the appropriate function at 1 and use that*/
	/* as the output */
	if (pB_name==1) {
		BUTO = polyeval(MTO_m,pB_name)
		BUUO = polyeval(MUO_m,pB_name)
		BUTE = polyeval(MTE_m,pB_name)

	} 
	else {
		BUTO = polyeval(polyinteg(w_BU_MTO,1),1) - polyeval(polyinteg(w_BU_MTO,1),pB_name)
		BUUO = polyeval(polyinteg(w_BU_MUO,1),1) - polyeval(polyinteg(w_BU_MUO,1),pB_name)
		BUTE = polyeval(polyinteg(w_BU_MTE,1),1) - polyeval(polyinteg(w_BU_MTE,1),pB_name)
	}
				
	st_numscalar("sc_BUTO", BUTO)
	st_numscalar("sc_BUUO", BUUO)
	st_numscalar("sc_BUTE", BUTE)	
								
	/* IT weights */

	/* Define the weights as 1 from 0 to pI*/
	/* There is no reference to 0 and pI here, but later when we take the */
	/* sum of the weights we take it from 0 to pI */
	IT_w = 1
		
	/* Integrate the weights from 0 to pI */
	/* This is equivalent to calculating the sum of the IT weights from */
	/* 0 to 1 */
	IT_sum = polyeval(polyinteg(IT_w,1),pI_name) - polyeval(polyinteg(IT_w,1), 0)
		
	/* Normalize the weights to sum to 1 by dividing the weights by the */
	/* sum from 0 to 1 */
	polydiv(IT_w, IT_sum, IT_w_new, temp)
		
	/* Multiply MTO, MUO, and MTE by the weights calculated above*/
	w_IT_MTO = polymult(MTO_m, IT_w_new)
	w_IT_MUO = polymult(MUO_m, IT_w_new)
	w_IT_MTE = polymult(MTE_m, IT_w_new)
		
	/* Integrate MTO, MUO, and MTE from 0 to pI */
	/* If pI = 0, then evaluate the appropriate function at 0 and use that*/
	/* as the output */
	if (pI_name==0) {
		ITTO = polyeval(MTO_m,pI_name)
		ITUO = polyeval(MUO_m,pI_name)
		ITTE = polyeval(MTE_m,pI_name)
	
	} 
	else {
		ITTO = polyeval(polyinteg(w_IT_MTO,1),pI_name) - polyeval(polyinteg(w_IT_MTO,1),0)
		ITUO = polyeval(polyinteg(w_IT_MUO,1),pI_name) - polyeval(polyinteg(w_IT_MUO,1),0)
		ITTE = polyeval(polyinteg(w_IT_MTE,1),pI_name) - polyeval(polyinteg(w_IT_MTE,1),0)
	}
				
	st_numscalar("sc_ITTO", ITTO)
	st_numscalar("sc_ITUO", ITUO)
	st_numscalar("sc_ITTE", ITTE)			
								
	/* IU weights */

	/* Define the weights as 1 from pI to 1*/
	/* There is no reference to 1 and pI here, but later when we take the */
	/* sum of the weights we take it from pI to 1 */
	IU_w = 1
		
	/* Integrate the weights from pI to 1 */
	/* This is equivalent to calculating the sum of the IU weights from */
	/* 0 to 1 */
	IU_sum = polyeval(polyinteg(IU_w,1),1) - polyeval(polyinteg(IU_w,1), pI_name)
		
	/* Normalize the weights to sum to 1 by dividing the weights by the */
	/* sum from 0 to 1 */
	polydiv(IU_w, IU_sum, IU_w_new, temp)
		
	/* Multiply MTO, MUO, and MTE by the weights calculated above*/
	w_IU_MTO = polymult(MTO_m, IU_w_new)
	w_IU_MUO = polymult(MUO_m, IU_w_new)
	w_IU_MTE = polymult(MTE_m, IU_w_new)
		
	/* Integrate MTO, MUO, and MTE from pI to 1*/
	/* If pI = 1, then evaluate the appropriate function at 1 and use */
	/* that as the output */
	if (pI_name==1) {
		IUTO = polyeval(MTO_m,pI_name)
		IUUO = polyeval(MUO_m,pI_name)
		IUTE = polyeval(MTE_m,pI_name)
	
	} 
	else {
		IUTO = polyeval(polyinteg(w_IU_MTO,1),1) - polyeval(polyinteg(w_IU_MTO,1),pI_name)
		IUUO = polyeval(polyinteg(w_IU_MUO,1),1) - polyeval(polyinteg(w_IU_MUO,1),pI_name)
		IUTE = polyeval(polyinteg(w_IU_MTE,1),1) - polyeval(polyinteg(w_IU_MTE,1),pI_name)
	}
				
	st_numscalar("sc_IUTO", IUTO)
	st_numscalar("sc_IUUO", IUUO)
	st_numscalar("sc_IUTE", IUTE)				
		
				
	/* LA weights */
		
	/* Define the weights as 1 from pB to pI*/
	/* There is no reference to pB and pI here, but later when we take the*/
	/*  sum of the weights we take it from pB to pI*/
	LA_w = 1
		
	/* Integrate the weights from pB to pI */
	/* This is equivalent to taking the sum of the LA weights from 0 to 1 */
	LA_sum = polyeval(polyinteg(LA_w,1),pI_name) - polyeval(polyinteg(LA_w,1), pB_name)
		
	/* Normalize the weights to sum to 1 by dividing by the sum*/
	polydiv(LA_w, LA_sum, LA_w_new, temp)
		
	/* Multiply MTO, MUO, and MTE by the weights */
	w_LA_MTO = polymult(MTO_m, LA_w_new)
	w_LA_MUO = polymult(MUO_m, LA_w_new)
	w_LA_MTE = polymult(MTE_m, LA_w_new)
		
	/* Integrate MTO, MUO, and MTE from pB to pI */
	/* If pB = pI, then evaluate the appropriate function at pB (or pI, */
	/* doesn't matter) and use that as the output */
	if (pI_name==pB_name) {
		LATO = polyeval(MTO_m,pI_name)
		LAUO = polyeval(MUO_m,pI_name)
		LATE = polyeval(MTE_m,pI_name)
	}
	else {
		LATO = polyeval(polyinteg(w_LA_MTO,1),pI_name) - polyeval(polyinteg(w_LA_MTO,1),pB_name)
		LAUO = polyeval(polyinteg(w_LA_MUO,1),pI_name) - polyeval(polyinteg(w_LA_MUO,1),pB_name)
		LATE = polyeval(polyinteg(w_LA_MTE,1),pI_name) - polyeval(polyinteg(w_LA_MTE,1),pB_name)
	}

	st_numscalar("sc_LATO", LATO)
	st_numscalar("sc_LAUO", LAUO)
	st_numscalar("sc_LATE", LATE)
			
	/* A weights */
		
	/* Define the weights as 1 */
	A_w = 1
		
	/* Integrate the weights from 0 to 1 - this is equivalent to taking */
	/* the sum of the weights from 0 to 1*/
	A_sum = polyeval(polyinteg(A_w,1),1) - polyeval(polyinteg(A_w,1), 0)
		
	/* Normalize the weights to sum to 1 by dividing by their sum */
	polydiv(A_w, A_sum, A_w_new, temp)
		
	/* Multiply MTO, MUO, and MTE by the weights */
	w_A_MTO = polymult(MTO_m, A_w_new)
	w_A_MUO = polymult(MUO_m, A_w_new)
	w_A_MTE = polymult(MTE_m, A_w_new)
		
	/* Integrate MTO, MUO, and MTE from 0 to 1 */
	ATO = polyeval(polyinteg(w_A_MTO,1),1) - polyeval(polyinteg(w_A_MTO,1),0)
	AUO = polyeval(polyinteg(w_A_MUO,1),1) - polyeval(polyinteg(w_A_MUO,1),0)
	ATE = polyeval(polyinteg(w_A_MTE,1),1) - polyeval(polyinteg(w_A_MTE,1),0)
	
	st_numscalar("sc_ATO", ATO)
	st_numscalar("sc_AUO", AUO)
	st_numscalar("sc_ATE", ATE)
		
	/* RIST weights */

	/* Define the unnormalized weights as all treated individuals */
	/* (baseline treated + treated compliers)*/
	/* We define the weights for the baseline treated and the treated */
	/* compliers separately because the treated compliers are obtained */
	/* using the LATE weights weighted by the fraction randomized in */
	
	/* This is the first piece of the RIST weights, for ALL the baseline */
	/* treated */
	RIST_w_1 = 1
		
	/* This is the second piece of the RIST weights, for the compliers, */
	/* weighted by the fraction that was randomized in */
	RIST_w_2 = s_pI_name
		
	/* Sum the weights for the baseline treated and compliers separately */
		
	/* Sum the weights for the baseline treated - from 0 to pB */
	RIST_w_1_sum = polyeval(polyinteg(RIST_w_1,1),pB_name) - polyeval(polyinteg(RIST_w_1,1), 0)
	
	/* Sum the weights for the compliers - from pB to pI, previously */
	/* weighted by the fraction that was randomized in and treated */
	RIST_w_2_sum = polyeval(polyinteg(RIST_w_2,1),pI_name) - polyeval(polyinteg(RIST_w_2,1), pB_name)
	
	/* Sum the weights for the baseline treated and the treated compliers */
	/* together to get the sum of the weights over the range from 0 to 1*/
	RIST_w_sum = RIST_w_1_sum + RIST_w_2_sum
		
	/* Normalize the weights for the baseline treated and the treated */
	/* compliers to sum to 1 alltogether by dividing by the sum of the */
	/* weights, but normalize them separately*/
	
	polydiv(RIST_w_1, RIST_w_sum, RIST_w_1_new, temp)
	polydiv(RIST_w_2, RIST_w_sum, RIST_w_2_new, temp)
		
	/* Multiply MTO, MUO, and MTE by the weights, separately for the */
	/* baseline treated and the treated compiers */
	
	/* MTO*/
	w_RIST_1_MTO = polymult(MTO_m, RIST_w_1_new)
	w_RIST_2_MTO = polymult(MTO_m, RIST_w_2_new)
	
	/* MUO*/
	w_RIST_1_MUO = polymult(MUO_m, RIST_w_1_new)
	w_RIST_2_MUO = polymult(MUO_m, RIST_w_2_new)
	
	/* MTE*/
	w_RIST_1_MTE = polymult(MTE_m, RIST_w_1_new)
	w_RIST_2_MTE = polymult(MTE_m, RIST_w_2_new)
		
	/* Integrate MTO, MUO, and MTE from 0 to pB for the baseline treated */
	/* and from pB to pI for the treated compliers. */
	/* We integrate the MTO, MUO, and MTE separately for the baseline */
	/* treated and the compliers, but then we add the result together to */
	/* get the outcome/effect for all treated individuals*/
		
	/*MTO*/
	RISTTO1 = polyeval(polyinteg(w_RIST_1_MTO,1),pB_name) - polyeval(polyinteg(w_RIST_1_MTO,1),0)
	RISTTO2 = polyeval(polyinteg(w_RIST_2_MTO,1),pI_name) - polyeval(polyinteg(w_RIST_2_MTO,1),pB_name)
	RISTTO = RISTTO1 + RISTTO2
	
	/*MUO*/
	RISTUO1 = polyeval(polyinteg(w_RIST_1_MUO,1),pB_name) - polyeval(polyinteg(w_RIST_1_MUO,1),0)
	RISTUO2 = polyeval(polyinteg(w_RIST_2_MUO,1),pI_name) - polyeval(polyinteg(w_RIST_2_MUO,1),pB_name)
	RISTUO = RISTUO1 + RISTUO2
	
	/*MTE*/
	RISTTE1 = polyeval(polyinteg(w_RIST_1_MTE,1),pB_name) - polyeval(polyinteg(w_RIST_1_MTE,1),0)
	RISTTE2 = polyeval(polyinteg(w_RIST_2_MTE,1),pI_name) - polyeval(polyinteg(w_RIST_2_MTE,1),pB_name)
	RISTTE = RISTTE1 + RISTTE2

	st_numscalar("sc_RISTTO", RISTTO)
	st_numscalar("sc_RISTUO", RISTUO)
	st_numscalar("sc_RISTTE", RISTTE)
				
	/* RISU weights */

	/* Define the unnormalized weights as all untreated individuals */
	/* (intervention untreated + untreated compliers) */
	/* We define the weights for the intervention untreated and the */
	/* untreated compliers separately because the untreated compliers are */
	/* obtained using LATE weights weighted by the fraction randomized out*/
	
	/* This is the first piece of the RISU weights, for the compliers, */
	/* weighted by the fraction that was randomized out */
	RISU_w_1 = s_pB_name
	
	/* This is the second piece of the RISU weights, for ALL of the*/
	/* intervention untreated */
	RISU_w_2 = 1
	
	/* Sum the weights for the intervention untreated and the compliers */
	/* separately */
		
	/* Sum the weights for the compliers - from pI to pB, previously */
	/* weighted by the fraction that was randomized out and untreated*/
	RISU_w_1_sum = polyeval(polyinteg(RISU_w_1,1),pI_name) - polyeval(polyinteg(RISU_w_1,1), pB_name)
	
	/* Sum the weights for the intervention untreated - from pI to 1 */
	RISU_w_2_sum = polyeval(polyinteg(RISU_w_2,1),1) - polyeval(polyinteg(RISU_w_2,1), pI_name)
	
	/* Sum the weights for the intervention untreated and the untreated */
	/* compliers together to get the sum of the weights over the full range from 0 to 1*/
	RISU_w_sum = RISU_w_1_sum + RISU_w_2_sum
					
	/* Normalize the weights for the intervention untreated and the */
	/* untreated compliers to sum to 1 alltogether by dividing by the sum */
	/* of the weights, but normalize them separately*/
	
	polydiv(RISU_w_1, RISU_w_sum, RISU_w_1_new, temp)
	polydiv(RISU_w_2, RISU_w_sum, RISU_w_2_new, temp)
	
	/* Multiply MTO, MUO, and MTE by the weights, separately for the */
	/* intervention untreated and the untreated compliers */
	
	/* MTO*/
	w_RISU_1_MTO = polymult(MTO_m, RISU_w_1_new)
	w_RISU_2_MTO = polymult(MTO_m, RISU_w_2_new)
	
	/* MUO*/
	w_RISU_1_MUO = polymult(MUO_m, RISU_w_1_new)
	w_RISU_2_MUO = polymult(MUO_m, RISU_w_2_new)
	
	/* MTE*/
	w_RISU_1_MTE = polymult(MTE_m, RISU_w_1_new)
	w_RISU_2_MTE = polymult(MTE_m, RISU_w_2_new)
	
		
	/* Integrate MTO, MUO, and MTE from pB to pI for the untreated */
	/* compliers and from pI to 1 for the intervention untreated. */
	/* We integrate the MTO, MUO, and MTE separately for the intervention */
	/* untreated and the compliers, but then we add the result together to */
	/* get the outcome/effect for all untreated individuals*/
		
	/*MTO*/
	RISUTO1 = polyeval(polyinteg(w_RISU_1_MTO,1),pI_name) - polyeval(polyinteg(w_RISU_1_MTO,1),pB_name)
	RISUTO2 = polyeval(polyinteg(w_RISU_2_MTO,1),1) - polyeval(polyinteg(w_RISU_2_MTO,1),pI_name)
	RISUTO = RISUTO1 + RISUTO2
	
	/*MUO*/
	RISUUO1 = polyeval(polyinteg(w_RISU_1_MUO,1),pI_name) - polyeval(polyinteg(w_RISU_1_MUO,1),pB_name)
	RISUUO2 = polyeval(polyinteg(w_RISU_2_MUO,1),1) - polyeval(polyinteg(w_RISU_2_MUO,1),pI_name)
	RISUUO = RISUUO1 + RISUUO2
	
	/*kp*/
	RISUTE1 = polyeval(polyinteg(w_RISU_1_MTE,1),pI_name) - polyeval(polyinteg(w_RISU_1_MTE,1),pB_name)
	RISUTE2 = polyeval(polyinteg(w_RISU_2_MTE,1),1) - polyeval(polyinteg(w_RISU_2_MTE,1),pI_name)
	RISUTE = RISUTE1 + RISUTE2

	st_numscalar("sc_RISUTO", RISUTO)
	st_numscalar("sc_RISUUO", RISUUO)
	st_numscalar("sc_RISUTE", RISUTE)	
						
}
end
